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Abstract 

Interpolation error estimates needed in common finite element applications using simplicial 
meshes typically impose restrictions on the both the smoothness of the interpolated functions 
and the shape of the simplices. While the simplest theory can be generalized to admit less 
smooth functions (e.g., functions in H^{il) rather than iJ^(r2)) and more general shapes (e.g., 
the maximum angle condition rather than the minimum angle condition), existing theory does 
not allow these extensions to be performed simultaneously. By localizing over a well-shaped 
auxiliary spatial partition, error estimates are established under minimal function smoothness 
and mesh regularity. This construction is especially important in two cases: L^(i7) estimates 
for data in W^'^{il,) hold for meshes without any restrictions on simplex shape, and W^'P{^1) 
estimates for data in W'^'P{fl) hold under a generalization of the maximum angle condition 
which requires p > 2 for standard Lagrange interpolation. 

Interpolation error estimates for the standard finite element spaces typically involve two re- 
quirements: conditions on the smoothness of the function being interpolated (data regularity) and 
constraints on the geometry of the elements of the mesh (shape regularity). Beyond the simplest 
interpolation error estimates (stated precisely in Proposition II. 2p . data or shape regularity can be 
substantially relaxed with some additional analysis. However the existing theory does not accept 
the weakest restrictions on both data and shape regularity simultaneously. 

Data Regularity and Average Interpolation Lagrange interpolation involves point-wise func- 
tion evaluation leading to data regularity restrictions associated with the requirements of the 
Sobolev embedding theorem. Alternatively, average interpolation can be used to handle less regular 
data |13 1 138 ^ 151 137 1 [TT]. Because average interpolation 'smears' local estimates over a neighborhood 
of simplices, error estimates rely more heavily on the geometry of the mesh, often needing meshes 
of bounded ply (defined in Section 11.31 see Figure [I]) so global summation can be performed. 

Shape Regularity and Angle Conditions For triangular meshes the classical error estimates 
are typically proved under a minimum angle restriction, but this condition is overly restrictive. 
Triangle geometry and interpolation error are actually related through the largest angle of the 
triangle [431 [6| [22] . (Convergence of the finite element method can occur for meshes with arbitrarily 
large angles [8l[T9], but the known counterexamples are finite element spaces containing subspaces 
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associated with a shape regular mesh.) The maximum angle condition can be generalized to higher- 
dimensional simplices but existing analysis gives stronger data regularity requirements [221 [Ml [39] . 
The linear tetrahedral element is especially important to this discussion. Shenk constructed a 
counterexample demonstrating that the expected Lagrange interpolation error estimates do not 
hold for a fairly innocuous family of narrow tetrahedra |39) : see Fig. [3{right). Analysis using an 
average interpolant eliminates this problem but again leads to new geometric restrictions on the 
ply of the mesh [Ij. 

Mesh Generation Geometric restrictions needed to prove interpolation error estimates are es- 
sentially the output requirements for a mesh generator. In 2D many mesh generation techniques 
aim to produce bounded aspect ratio triangulations (although a notable exception f30] only removes 
large angles), but aspect ratio guarantees are impossible when the input (i.e., the boundary of the 
domain) contains small angles [30]. To admit arbitrary input several mesh generation strategies 
have been developed which allow a few poor quality triangles near small input angles |29[ l33 [[30ll36] . 
Guided by the interpolation theory the most commonly used strategy p9] provides provable upper 
bounds on the largest angle and the ply of the mesh (as long as the input does not include high 
degree vertices). This effective approach does not extend to 3D and guaranteed mesh generation 
algorithms do not yield a bound on the ply [251 1121 HIl Hffl IMl [36]. Figure [Deleft) demonstrates 
the high ply construction in the analogous 2D algorithm. Failures of these algorithms are difficult 
to observe in practice since the mesh ply only grows logarithmically in the (inverse of the) size of 
the smallest angle, but they represent a theoretical disconnect between output guarantees of the 
mesh generator and the input requirements of numerical methods. Beyond provably-correct mesh 
generation algorithms, simplified mesh generation algorithms are often practically successful, and 
allowing high ply meshes extends the range of inputs producing acceptable results. Uniform refine- 
ment of a structured grid in polar coordinates is perhaps the simplest situation yielding a high ply 
node; see Figure [2][right). Providing the loosest set of geometric requirements in the interpolation 
theory gives mesh generators the most flexibility in practical applications as well as new challenges 
to improve existing guarantees under weaker requirements. 

Outline This paper demonstrates that there is no inherent trade-off between data regularity and 
shape regularity in interpolation of less regular Sobolev functions on simplicial meshes by developing 
a theory accepting minimal data and shape regularity requirements. The key insight is that error 
estimates cannot be localized over the meshes of high ply; rather a simpler auxiliary spatial partition 
must be used. Before describing this construction. Section [Tj contains the necessary preliminary 
discussion including a more formal exposition of the prior results. Error estimates for the average 
interpolant under uniform size and quality meshes are shown in Section [2l A generalization to 
non-uniform meshes is given in Section [3] before some concluding remarks in Section [H 

1 Preliminaries 

The statement x < y is used to mean there exists a positive constant C such that x < C y when the 
details of the constant are unimportant. For nearly all estimates discussed, this implied constant 
will be independent of simplex/mesh shape regularity; whenever this is not the case, it will be 
explicitly noted. 
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Figure 1: Two refinement strategics that do not produce bounded aspect ratio meshes, (left) A 
refinement strategy creating a mesh with bounded ply and maximum angle: all vertices have degree 
six and all angles are acute, (right) A refinement strategy creating a high ply mesh without large 
angles: all angles are smaller than 135 degrees: in the limiting case, the triangles in the corners 
contain angles very near 135°, 45° and 0°. 





Figure 2: (left) A quality mesh generated for an in- 
put with a very small angle. Certain mesh generation 
strategies can produce high ply meshes at vertices (or 
along edges in 3D) near small input angles, (above) 
"Uniform" refinement in polar coordinates leads to 
a high ply at the origin. 
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1.1 Balls, Cubes, and Simplices 

The ball of radius r centered at x G M"^ is the set of points with Euchdean distance from x less 
than r: 

B(x,r) := {y : |x - y| < r}. 
A d-cube or cube is a closed set Q C M'^ of the form 



Q 



h 



h 



h h 
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The point a = (ai, . . . , a^) is called the center of Q and h is called the si^e of Q, denoted size((5)- 
A simplex is the convex hull of d + 1 non-coplanar points (called vertices), {vjj^^Q: 



K 



^tiVi : ti > 0,^ti = 1 



. i=0 



1=0 



We let /li^ denote the length of the longest side of K and pK denote the radius of the largest 
sphere inscribed in K. Following Jamet's shape regularity condition [22], the angle 6k is defined 
to characterize the quality of a simplex K by the formula, 



Ok '■= maxmin arccos^ • Uj 
1^1=1 i 



(1) 



where {uj} is the set of unit vectors parallel to the edges of the simplex K; see Fig. [3l This 
angle measures how far from coplanar the set of edges is, so we call 9k the coplanarity measure. 
The subscript on h, p, and 9 is often omitted when the simplex in question is apparent. In two 
dimensions 9k is half the largest angle of the triangle, and thus the set of triangles with coplanarity 
bounded away from tt/2 satisfy the maximum angle condition in the sense that all angles are 
bounded away from vr. In three dimensions, a more intuitive generalization of the maximum angle 
condition defined by defined by Kfizek [21] requires that both face and dihedral angles of the 
tetrahedron are bounded away from vr. This is also equivalent to bounding coplanarity away from 
7r/2 as shown in the proposition below. 

Proposition 1.1. LetipK denote the maximum of all face and dihedral angles of (three-dimensional) 
tetrahedron K. For tetrahedra, bounds on coplanarity and face/dihedral angles are equivalent, i.e., 
if 9k is bounded away from 7r/2 then ipK is bounded away from vr and vice versa. 

The proof, given in Appendix |Al relies on a characterization of tetrahedra satisfying the maxi- 
mum angle condition given in [2]. 



1.2 Sobolev Spaces 

We will distinguish between the notation ^ and di. The former will be used to denote differen- 
tiation with respect to a particular variable while the latter represents differentiation with respect 
to a function argument. For example, the following notation is used for an application of the chain 
rule: 

— [u(2x + 3y)] = 2diu{2^ + 3y). 
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Figure 3: Coplanarity 9 for for several simplices with vector ^max which maximizes the expression 
in ([T]) labeled, (left) In 2D the coplanarity is half the largest angle, (center) The coplanarity of a 
nearly flat "sliver" tetrahedron is nearly 7r/2. (right) The simplices used in a counter-example due 
to Shenk |39j have coplanarity bounded away from 7r/2 even as the edge length e approaches 0. 

The Sobolev semi-norms and norms on functions over an open set VL are defined by 



where the first summation is taken over multi-indices j3 and differentiation with respect to multi- 
indices has the standard definition, d^u = d^^d^^ ...u. The VF'^'^-norm is simply the standard 
L^'-norm. 

Throughout this paper, C M*^ is a bounded set and an extension domain for any Sobolev 
spaces used; i.e. there is a continuous, linear operator Eq, : VF'"'^'(0) — )• W"^'P{M.'^) such that 
u{x) = Equ{x) for all x € fi; see [31 p. 83] for technical results concerning extension domains. 
Without ambiguity, the functions u and E^^u will not be distinguished: for x ^ fi, u{x) denotes 
Efiu{x). Likewise, the estimate |ii|i4/m,p(igd-) < l''^lvK'"'P(r2) used freely. 

1.3 Lagrange Interpolation on Simplicial Meshes 

T always denotes a simplicial mesh of The ply of the mesh T is the maximum number of mesh 
simplices intersecting a single simplex. Ply restrictions are often implicit in interpolation error 
estimates since a mesh with bounded aspect ratio has a bounded ply. (In fact, the weaker restriction 
of bounded circumradius to shortest edge ratio is sufficient to ensure bounded ply [441 [2T].) 
Let Ilfc denote the fc-th degree Lagrange interpolant on a simplicial mesh. 



where 4> are the piecewise nodal basis functions on T- While this paper only considers linear 
interpolants {k = 1), the higher order cases are included in the statement of classical estimates for 
smooth functions. 

1.4 Prior Interpolation Error Estimates 

The Bramble-Hilbert Lemma [71 along with scaling properties of Sobolev semi-norms leads 
to the following classical error estimate for the Lagrange interpolant; see [121 \TE\ [9] for complete 
details. 




and 




0<fc<m 
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Proposition 1.2. If k + 1 > ^ and < m < k, then for all simplices K and all u G W''^^'P{K), 

fl fc+1 

For simplices with bounded aspect ratio, h and p are proportional and thus the factor 
reduces to h^~^^~'^. The estimate in the L^-norm (i.e., the m = case) is particularly nice: the 
term in the denominator disappears and the estimate is independent of the shape of the simplex. 

The first generalization involves weakening the data regularity requirement, i.e., the restriction 
u S VF^'^(r2). A function u E W'^'^{^) (but not in VF^'^(O)) is not necessarily continuous so 
an alternative to Lagrange interpolation must be considered. Clement first described such an 
interpolant and proved an estimate which corresponds to the (A; = 0)-case of ([2]) [13]. Scott and 
Zhang proved a similar estimate based on an average interpolant [38]: specifically interpolant values 
at mesh vertices are defined to be averages over an arbitrarily selected adjacent simplex. A distinct 
advantage of the Scott and Zhang interpolant is that coefficients of basis functions corresponding to 
boundary nodes only depend upon the boundary data. An alternative approach based on averaging 
over balls has been developed which extends to vector-valued functions in the discrete de Rham 
complex [5| I371[TT] : see Figure[4]for a comparison of the averaging regions. The resulting interpolant 
of Christiansen and Winther also matches homogeneous boundary data selecting an appropriate 
extension of the function outside the domain. In the simplest case, each of these constructions 
reduces to the following theorem. 

Theorem 1.3 (Specialized from [HUlIl], etc.). Suppose T has a bounded aspect ratio and let Uq 
be an average interpolant. Then for all u G T^^'^(O), 

(3) 

where K is the union of K and its neighboring simplices. 

Unlike the (m = 0, /c = l)-case of Proposition II. 2^ Theorem 11.31 requires a bounded aspect ratio 
mesh ensuring two things. First, K and its neighboring simplices have roughly the same size so 
that diam(/C) w diam(K). Second, K has a bounded ply so these local estimates ([3]) over each 
simplex can be summed to a global estimate. Our aim is to remove these restrictions on simplex 
shape to match the m = case of the classical estimate ([2]) which is independent of simplex shape. 

The second generalization of Proposition 11.21 weakens the shape regularity requirements. For 
this purpose, Jamet developed an improved estimate for sufficiently regular data. 

Theorem 1.4 (Jamet, [22]). Ifk + l-m>^, then for all u G W^+'^^PiK), 

The restriction k + 1 — m > ^ prevents Jamet 's result from applying in certain cases including 
(perhaps the most important cases of) linear interpolation on triangles and tetrahedra: k = m = 1 
when d = 2 or d = 3. The analysis of Babuska and Aziz [6] yields the result for the linear 
triangular element (although care must be taken to derive the correct geometric factor l/cos0 
using this approach |23[ [T8\ [28l I34j). Shenk generalized the approach of Babuska and Aziz to 
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higher dimensions and, more importantly, demonstrated that @ does not hold for the linear 
Lagrange interpolant on tetrahedra [39]. Shenk's example involves particularly simple geometry: a 
family of tetrahedra depicted in Figure [3] with vertices (0, 0, 0), (1, 0, 0), (0, 1,0), and (0, 0, e) where 
e approaches zero. Demonstrating the error estimate for an analogous class of triangles (i.e., 2D 
simplices) is the crux of the argument of Babuska and Aziz. 

An analysis by Acosta demonstrated that expected error estimates can be extended to certain 
narrow tetrahedra (including those in Shenk's example) by using average (rather than Lagrange) 
interpolation [1|. This method does allow data and shape regularity to be weakened simultaneously, 
but the analysis is inherently anisotropic: locally it must be possible to affinely transform the mesh 
into a shape-regular mesh. (This is a standard requirement in an anisotropic framework; e.g., 
[31[I71[20].) The resulting estimates expect the mesh to have a bounded ply so that local estimates 
can be summed to form global error estimates. Note that the meshes in Figure [2] do not satisfy the 
requirements of [I]. 

We will develop interpolation error estimates based on average interpolation under weakened 
data and shape regularity requirements. Emphasis will be placed on avoiding any unnecessary 
element shape requirements which, most notably, allows meshes with unbounded ply (as is the 
case for Lagrange interpolants in Theorem ll.4p . The result will be an interpolation theory which 
extends Jamet's minimal shape regularity requirements to less smooth data including the linear 
tetrahedral element. 

2 Uniform Error Estimates 

To avoid the limitations associated with pointwise function evaluation, an average interpolant is 
defined by mollifying the data before Lagrange interpolation: 

Ehu{x) := UiMhu{x). (5) 

In the next subsection, the standard mollification procedure A^h is defined and a number of needed 
properties are given. 

Remark 2. A. This interpolant is used in [5} 137^ ITT], but in that context, the mollification radius 
was selected to be small enough so that the averaging region stays within the neighboring triangles. 
Since we do not assume aspect ratio bounds, we cannot use this property. Acosta uses essentially 
an anisotropic variant of this construction [Tj; see Figured! 

2.1 Mollification and Integration 

Mollification is a common tool in analysis for producing smooth approximations of functions in 
spaces. After defining the procedure, a few common facts about mollification are listed. The proofs 
of these and similar results can be found in a number of textbooks, but are given in the appendix for 
comparison to the later sections. For some examples of these kinds of technical details, especially 
those for more regular data and including explicit dependence on on the mollification radius, see 
[all p. 191], [261 P- 58], [271 p. 98], and [16, p . 201]. 

Let p G C~(M'^) such that (i) supp(p) C B{0,1), (ii) / p(x) dx = 1, and (iii) p(x) > for ah 
x. Let ph{^) = 7^p(f )• Then, (i) supp(p/,) C B{0,h), (ii) /y9/j(x)dx = 1, (iii) p/,(x) > for aU x, 

and (iv) ||Vpft||^i^]gd-) < ||Vp||j;^i^]gd). The most commonly used mollifier is p(x) = ce^-'^/^^^H^H ^ 
where the constant c is chosen so (ii) holds. We also require that p is radially symmetric. While 
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Figure 4: Averaging regions for the interpolants of (left to right) Clement, Scott and Zhang, 
Christiansen and Winther, Acosta, and this paper. The first three were developed for shape reg- 
ular meshes. The final two are shown on an anisotropic mesh to highlight the differences in the 
construction. While the Acosta interpolant requires well-defined anisotropy in the mesh, our inter- 
polant makes minimal assumptions on the mesh quality. 

symmetry is not usually required, the common mollifier is symmetric and we explicitly use the 
property in one step of the upcoming arguments. 

A function is mollified by convolving it with p^: given u G L^(W^), define 

Mhu{x.) := / ti(x - y)pft(y) dy = / n(x - /iz)p(z) dz. 

The two representations above are equivalent, resulting from a simple change of variables, y = hz, 
but the supports of the integrands are different. If a function is smooth enough (specifically, it 
belongs to W^'^{Q) in the proposition below) then the difference between the function and its 
mollification can be bounded in terms of the mollification radius. 

Proposition 2.1. If p > 1, then for all u G W^'P{Cl), 

\\U- Mhu\\LP(^^d) < • 

The next proposition demonstrates how the smoothness of the mollifier depends on the molli- 
fication radius through a bound on the L^-norm of a mollified function in terms of the L'^-norm of 
the original function. 

Proposition 2.2. Ifl<q<p<oo and Q is a hounded set, then for all u G L'^{Q), 

I \^hu\\LP(Q) <h~\ |m| Il,(q^) (6) 

where Qh ■= UxeQ -^(^' 

Remark 2.B. Since i-Mhu) = -M-h Proposition 12.21 also applies for higher Sobolev 

norms; specifically, ([6]) implies that 

^ q — p 

\MhU\wk,P(^Q)<h vq \u\^k,^Q^y (7) 

Finally a standard estimate of the L^-norm of a function in terms of the L'^-norm is stated for 
q larger than p. This can be proved with either Jensen's or Holder's inequality. 

Proposition 2.3. Ifl<p<q<oo and Q is a bounded set, then for all u G L''{Q), 

\\u\\lp(^q) < \Q\^ M\li{Q) ■ 
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Figure 5: For an example domain, the cubes in Q are shaded. 
2.2 Error Estimate 

Next we turn to removing the integrabihty condition in Jamet's estimate dH). Again for simphcity 
we first consider meshes in which simplex size and quality (measured by coplanarity, not aspect 
ratio) are uniformly bounded and then relax these assumptions in Section [3l 

Theorem 2.4. Suppose T contain no edges longer than h and let 9 be the maximum simplex 
coplanarity. Then for k G {0, 1} and for all u E W''~^^'P{i^), 

Proof. First the triangle inequality is applied using the intermediate function AihU, 

\u - Ehu\^,k,p(^Q) <\u- Mhu\y^,k,p(^Q) + \MhU - Ehu\^,k,p(^Q) ■ (9) 

The first term of Q is independent of the triangulation T and is estimated using Proposition 12.11 

\u - Mhu\^k,p(^^) < h |n|^^fc+i,p(j^) . (10) 

Next the second term of ([9]) is addressed: 

\MhU - Ehu\^,k,p(^^) = \MhU - UiMhu\y^,k,p(^^) ■ (11) 

The estimate is established locally, but not over the individual simplices of T. This is important 
since averaging 'smears' the estimate on one simplex over its neighbors. Since we allow for high 
ply meshes, it is not possible to sum these local estimates to yield the global result dH). Rather, 
the estimate is localized over a simple artificial (and bounded-ply) decomposition, a covering of Q 
by cubes. 

In a lattice of cubes of side-length /i, let Q = {Qi} be the subset which intersect fi; see FigureEl 
The term pT]) is estimated on a single cube Qi G Q. Let Ri and Si denote the cubes of side-length 
3/i and 5/i, respectively, centered around Qi, and let ICi be the set of simplices in the triangulation 
Th which intersect Qi] see Figure [H 
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Next a value q > d is selected arbitrarily with the intention of applying Theorem 11.41 (since we 
are in the setting A; + l — m = l). Applying Proposition 12.31 followed by the definition of /Cj yields, 

\MhU - UiMhu\l,,^^^Q^^ < \Q^\^ \MhU - UiMhu\l,k^^^Q,^ 



Theorem 11.41 can now be applied to each term of the summation on the right hand side and 
Proposition 12.21 is used to return the estimate to the correct norm: 



{cos9y 



vy*+i.p(5i) 



Since \Qi\ = h'^, we conclude that 



\MhU - UiMhu\^yk,p,Q. < ^ . • (12) 

(cos a) 



Finally these estimates can be summed to get the global estimate: 

\MhU - UiMhuf^,^^^^^ < J2 \^hu - UiMhuf^.^^^Q^^ 



~l^(cos0)'=) ((cos0)'=) ''^''^'=+^0) • (13) 

The final inequality holds since each cube Qi belongs to at most 5'^ cubes Sj. □ 

Remark 2.C. The essential difference between this analysis and that of Acosta yy is the use of 
the artificial (and bounded-ply) mesh Q for establishing the local estimate (jl2p . 

Remark 2.D. Looking closely at the inequality ()13p reveals that the extension domain property 

has been used to estimate |^^lvi/2>p(U5i) I^lvy2.p(rj)- 

Remark 2.E. In the W case of ([8]), i.e, fc = 0, the estimate does not depend on simplex shape 
which agrees with the classical result for more regular data. 



3 Locally Quasi-Uniform Error Estimates 

Extending Theorem 12.41 to admit non-uniform meshes requires a precise characterization of the 
mesh size h at each point in the domain. We emphasize that the size function does not necessarily 
reflect the size of the triangles in the mesh. Rather it is only an upper bound. Because the average 
interpolant depends on a local neighborhood, the sizing function should not vary too rapidly so 
that there is roughly a unique size (up to a constant factor) associated with the neighborhood. 
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Figure 6: Diagram for the proof of Theorem 12.41 Triangles belonging to /Cj are shaded. 



Before proving error estimates, we outline precise restrictions on these sizing functions. Then 
based on the sizing function, variable-radius mollification and compatible cubic coverings are de- 
fined. As before, the final interpolation operator is a composition of the (modified) mollification 
operation with the usual Lagrange interpolation operator Hi: 



3.1 Sizing Function 

Using Lagrange interpolation for smooth (enough) functions, mesh size is captured locally by the 
length of the longest edge of a given triangle. For average interpolation, function values over 
neighboring triangles must be considered and thus the size of the averaging region also plays a 
role in the error estimate. On bounded aspect ratio meshes this is not an issue since the longest 
edges of any pair of neighboring triangles are within a constant factor [31] . To analyze average 
interpolants for meshes without aspect ratio bounds, we decouple the mesh sizing function from the 
actual size of the triangles. Rather, the sizing function will be a smooth (enough) function which 
provides an upper bound on the size of nearby triangles in the mesh. The result will be a sizing 
function which can be used in the error estimate despite admitting meshes with adjacent triangles 
with dramatically different sizes. 

Letting Csamg be a positive constant, we consider the class of sizing functions /i : ^ (0, oo 
with bounded first and second derivatives: 



S/iit(x) := UiMhu{x.). 



(14) 




(15) 




(16) 



We call a triangulation T of Q compatible with sizing function h if for all x E fi, /i(x) > /ij'(x) where 
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Figure 7: Comparison of two (bounded aspect ratio) meshes associated with sizing functions /i(x) = 
^ |x| + 0.01 (left) and /i(x) = |x|^ + 0.01 (right) over the square [—1,1]^. In each case, both the 
full mesh and the center region zoomed four times are shown. These meshes demonstrate the 
differences between conditions (jl5p and (jl6p : meshes only satisfying the former condition can 
grade more quickly away from zero. Note: the meshes do not satisfy (jlSp with the particular choice 
of constant, but this exaggerated example is given to more clearly display the effect. 

r(x) denotes the triangle containing x0 In this setting, we will establish uniform interpolation 
estimates for any sizing function with a compatible mesh as long as the sizing function satisfies 
(I15p and (I16p and the mesh satisfies the maximum angle condition (in the sense of coplanarity) . 

Remark 3. A. The particular choice of constant in (I15p is not essential. However, this value is 
very convenient for localizing future estimates to one layer of neighboring cubes in a 2:1 balanced 
cubic partition described in Subsection 13.31 Selecting a larger constant in (fT5]) will require analysis 
of more levels of the cubic partition. 

Remark 3.B. Restriction (I15p is implicitly present (although the constant is often weaker) in 
the context of shape regular meshes. For a shape-regular mesh T and letting P{T) be the union of 
triangles in T which intersect triangle T, we define 

hi'Vi) = max hx, 
TeT,Vi£P{T) 

This definition of h can be extended to using linear interpolation over the triangulation and it 
immediately follows that the triangulation T is compatible with h. Moreover, properties of well- 
spaced point sets ensure that \Vh\ can be bounded depending only on the mesh shape regularity [STj 
144} I21j . In other words, a shape regular mesh can be viewed as a mesh that is compatible with a 
Lipschitz sizing function, and thus (jl6p is the only extra condition on the sizing function we require 
to handle less shape-regular meshes. 

Remark 3. C. Restrictions (llSp and (I16p do allow for non-uniform meshes. The sequence of sizing 
functions /in(x) = |x|^ -|- can be paired with a sequence of meshes with progressively smaller 
elements near the origin. However, (|16p does pose some additional limitations on how quickly the 
sizing function can approach zero: in the shape regular case common adaptive refinement matches 
a sizing function /i„(x) = |x| + 2^ or /in(x) = min ||x| , which is not smooth enough for (fT6]) . 
See Figure [71 

Remark 3.D. Restrictions (llSp and (I16p are somewhat weaker than those used in [Tj where there 
are restrictions essentially of the form /i(x) < V/i(x). In that setting Gronwall's inequality prevents 



^Without causing confusion this is a slight abuse of notation using h{x) to denote the sizing function and hx to 
denote the length of the longest edge of triangle. This is done to emphasize the fact that the sizing function plays 
the role of the mesh size in the upcoming estimates. 
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such sizing functions from being truly non-uniform. Using the techniques of this paper, it should 
be possible to recover Acosta's anisotropic results in a similar relaxed setting. 



3.2 Variable-Radius Mollification 

We now consider the generalization of mollification by allowing the mollification radius to vary over 
the domain. Given u G L^{W^) and sizing function h, define 

Mhu{x.) := / u(x - y)ph(^) (y) dy = / n(x + /i(x)y)p(y) dy. 

The only difference between this definition and the standard mollifier is that the radius /i(x) 
is spatially varying. We consider an isotropic mollification region (and thus produce isotropic 
error estimates). The prior construction of Acosta uses an anisotropic mollifier, i.e., the term 
(/ii(x)yi, /i2(x)y2, /'-3(x)y3) replaces h{x)y in the above definition [1]. First we note a simple prop- 
erty of the variable radius mollifier. 

Proposition 3.1. Let p > 1 and Q C M.'^, and define Qh '■= IJ^^q -B(x, /i(x)). Then for all the 
sizing functions satisfying [T5\) and for all u G L'P{Qh), 

I |n(x + /i(x)y)|Pdx<||n||P,(Q^). (17) 

Proof. This results from a simple change of variables, z = x + /i(x)y. The Jacobian of this 
transformation can be bounded by y-^^ . By assumption ([15]) and since y G B{0,1), 1 -|- 

V/i(x) • y G [1 - l/(16\/d), 1 + l/(16\/d)] producing 1^7^). □ 

Extensions of Propositions 12.11 and 12.21 are necessary to prove the non-uniform error estimates. 
Since Remark 2.B no longer applies, i.e.. 



we prove these lemmas directly in the correct Sobolev norms. 

Lemma 3.2. Let p > 1 and Q he a cube. Then for k G {0, 1} and for all u G W'^'P{Q), 

\u - < i^Q ||u||^yfc+i,p(Q^) (18) 

where Hq = maxxgg ^(x), and Qh = UxeQ ^(^' ^(^))- 

Proof. The A; = case of inequality (jlSp follows from a very similar argument to the proof of 
Proposition 12. 1[ Letting w ^ L^ (M^) be an arbitrary function, 

{u{x)—Mhu{x))w{x)\ dx 

= ^ (^1^^^ ^(^(x - sh{x)y)) ds^ p{y) dyj w{x) 

< |/i(x)y • Vu(x - s/i(x)y)p(y)'u;(x)| dxdy ds. 

Jo JIRd Jq 



dx 
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Next we note that |/i(x)y| < Hq and then bound the inner integral independent of y by expanding 
the set Q: 

\{u{yi)-Mhu{^))w{-x.)\ dx 

Hq i I p(y) / |Vn(x - s/i(x)y)?x;(x)| dxdy ds 

JO JR'i Jo 



< 

< Ho \u\ 



Note: the construction J p(y) dy = 1 was used to ehminate the integral in the y variable. Selecting 
w = (u(x) — A^/jm(x))^ ^ completes the result. 

Next, the k = 1 case of inequality p8|) is addressed. We begin by computing -£^A4hu{x) 



directly: 



_d_ 

dxi 



Mhu{x) = Mhdiu{x) + / p{y)dih{x)Vu{x + /i(x)y) • y dy. 



Then letting g{x.) := J^a p{y)dih{x)Vu{x + /i(x)y) • y dy, 
^ {u-Mhu) 



dxi 



LP{Q) 



< \ \diu- Mh{diu)\\^^^Q^ + ||5|Ilp(q) 



The first term can be estimated using the argument applied to (jlSp yielding, 



(19) 



(20) 



(21) 



It only remains to estimate ||9'||2,p(q)- Let R^|_ denote the halfspace with positive first coordinate, 
Xi > 0. First g is rewritten in terms of second derivatives of u: 

g{x) = dih{x) I p(y) Vn(x + /i(x)y) • y dy 

= dihix) [ p{y) (Vn(x + /i(x)y) - Vu{x - h{x)y)) ■ y dy 

= ^^h{x) j^^ p{y) J ^ ^ Vn(x + s/i(x)y) ds^ ■ y dy. 

Our selection of a radially symmetric mollifier is important above as the property p{y) = p{—y) 
has been applied. 



\g{x)w{x)\ dx< / / [' |a,/i(x 
Jo Jr+ J-1 



)|p(y) 



d 

—\/u{x. + s/i(x)y) 

OS 



\y\ \w{x)\ dsdydx 



< / \s\ 



p{y) / |/i(x)| |D2^(x+s/i(x)y)| |u;(x)| dxdyds, 

Jq 



where D'^u denotes the Hessian of u, the matrix of second derivatives. Next, Holder's inequality 
can be applied: 



J l5(x)w(x)| dx < Hq \u\w2,j,(^q^) \H\lp'(q) ■ (22) 
Finally setting w{'k) = g((x)P~^, and then combining (|22p with (j2ip completes the result. □ 
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The next step is to develop an extension of Proposition 12.21 for the variable-radius mollifier. 



Lemma 3.3. Ifl<q<p<oo and Q is a bounded set, then for all the sizing functions satisfying 
and and for all u G W'^-'i{Q), k G {0, 1, 2}, 



\^huyk,p(^Q) < hg \Qh\ " IQl" \\u\\wk.,(Q^) (23) 

where hq = minxeg /i(x) and Qh ■= UxeQ -^(^' H^))- 
Proof. We begin with the case A; = 0. For all x & Q, 

\Mhu{^)\ < ||/?||icx,md) / |u(x + /i(x)y)| dy 

< ll/0|lrcx)md^ / /l(x)~'^ |u(z)| dz. 
JQi, 



The previous inequality results from a change of variables z = x + /i(x)y. (Note, x is fixed in this 
transformation.) Holder's inequality gives, 

\Mhu{^)\ = IIHlL-(Rd) h-'^lQhl^/"' \\u\\mQ^^ . (24) 

Integrating (PH) in x over Q produces in the A: = case. 

Next, we turn to the case k = 2; the k = 1 case follows from very similar arguments. 



d 



2 



dxidxj 



Mhu{:ic) = / p(y) 



didju{x + /i(x)y) + 

9i/i(x)V {dju{x + /i(x)y)) • y + 

djh{x)V {diu{x + /i(x)y)) • y + (25) 
dih{x)djh{x)y ■ D\ (x + /i(x)y)) y + 



didjh{x)V (u(x + /i(x)y)) • y 



dy 



The first term is identical to the k = case using the function didju in the place of u. In the 
second, third, and fourth terms, the same argument holds after we recall (jl5p and the fact that 
|y| < 1. While the argument is the same, the final term yields a few slight differences: (I16p is 
required and the resulting estimate is in terms of Klvi/i'P(Qh) ^^^her than 1^1 vi^2,p(q^) seen in all the 
other cases. This causes the right-hand side of ()23p to contain the full Sobolev norm, rather than 
the semi-norm. □ 

Remark 3.E. In Lemma 13.31 ()16p is only necessary for the case k = 2. 

Remark 3.F. Lemmas 13.21 and 13.31 will be applied for cubes Q such that hq w Hq size{Q) 
size(Q/i). In this setting the estimates yield the same factors of '/i' as seen in the analog with fixed 
mollification radius. 
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Algorithm 1 Cover 17 by cubes proportional to h. 



Require: Q = {Q} where Q is a cube such that size(Qj) > 4maxxen 
while 3Qi G Q such that size{Qi) > Smaxxeg^ h{x) do 
Remove Qi from Q 

Subdivide Qi into 2^ new cubes of size 
Insert all of the new cubes which intersect Q into Q. 
end while 



Figure 8: An example cubic covering of a domain for some non- uniform sizing function. 
3.3 Cubic Covering of 

An important idea in Section [2] was the use of a cubic grid of size h to cover the domain. In the 
non-uniform case we seek to cover with cubes such that each cube has side length proportional 
to /i(x) for X G Q; see Figure [8] for an example. This is possible due to restriction ([15]) and can be 
achieved through the simple algorithm given in Algorithm [TJ 

First we assert that any neighboring cubes produced have proportional sizes. 

Proposition 3.4. Let Q = {Qi} be the set of cubes resulting from Algorithm{ll If Qi n Qj ^ 0, 
then size{Qi) < 2size{Qj). 

Proof. Suppose Qi and Qj are adjacent cubes such that size{Qi) < Let Xj G n Qj and 

let Xj G Qj. Then, 

Thus Qj is eligible to be split by Algorithm [T] and we conclude that upon termination of the 
algorithm the lemma holds. □ 

The next proposition shows that for a compatible triangulation the simplices which intersect a 
particular cube are contained well inside the cube and its neighboring cubes. 

Proposition 3.5. Let Q be the set of cubes resulting from Algorithm [IJ Then for all x £ Q, 
4/i(x) < size(Q) < 16/i(x). 
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Proof. Let Q be any cube resulting from Algorithm [T] and let Q be the cube which was subdivided 
to form Q. Then 

2size((5) = size((5) > 8max/i(x) > 8max/i(x). (26) 

There exists x € Q such that /i(x) > '^'^"g'^^ since Q is not split further by the algorithm. By 
requirement p^ . for any y £ Q, 

Combining (j26p and (j27p yields the result. □ 

Lemma 3.6. Let triangulation T be compatible with h and let Q = {Qi} be the set of cubes resulting 
from AlgorithmUl Then for all Q £ Q, 

U U i?(x,Mx))c5.:= U Q„ 
where K,{Q) denotes the set of simplices intersecting Q. 

Proof. Let K G IC{Q) and let x £ K. Then x belongs to Q or a cube neighboring Q since by 
Propositions 13.41 and 13.51 dist(x, Q) < Hk < size((5)/4 < size{Q)/2 for any neighboring cube Q. 
Let Q denote the cube containing x and let y G -B(x, /i(x)). Also let z denote the nearest point in 
Q to X. The result is shown in two cases depending on the size of Q. 
Case 1 : size{Q) = size{Q) or size(Q) = size((5)/2. Then 

, , , , size(Q) size(Q) size(Q) 
|y - z| < |y - x| + |x - z| < /i(x) + Hk < + = 

By the balance property (Proposition 13. 4p . this implies that y G Ug^nQ^eQi- 
Case 2 : size{Q) = 2size(Q). Then 

|y - z| < |y - x| + |x - z| < /i(x) + hK 

^ , , , |z — x| size(Q) ^ . /I 1 

< h{z) + J jJ- + ^ < size(Q) - + 



wVd 4 ' ' V2 64\/d, 

The balance property (Proposition 13. 4P ensures that y must belong to a cube of size greater than 
or equal to size{Q). Since |y — z| < size{Q), y G Ug^nQT^fl Qj- '-' 

3.4 Error Estimate 

Now we have the necessary tools to prove an analog of Theorem [23] for a non-uniform sizing function 
and compatible triangulation. 

Theorem 3.7. Let k G {0,1}. For sizing functions h satisfying U5\} and il6\} paired with any 
compatible mesh T, let Q be the set of cubes produced by Algorithmic Then for all u G VF^'^(r2), 

where hi = maxxeg^ h{x) and 9* = maxrpnSi^ii ^T- 
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Note that this result is uniform over the class of sizing functions satisfying (jlSp and (jl6p . i.e., 
the inequality depends only upon Cgizing and not the specific sizing function or mesh selected. 

Proof. The proof structure is very similar to Theorem 12. 4i First the estimate is divided using the 
triangle inequality and the intermediate function 'EhU. 

\u - H,,u|^^fc,p(Q^) < |u - Mhu\y^,k.p(^Q^^ + \MhU - Ehu\^k^^^Q^^ . (29) 

Let Si denote the union of Qi with its neighboring cubes. By Lemmas 13.21 and 13.61 

Next the second term in ()29p is estimated. Select a value q > d. Proposition 12.31 and Theorem 11.41 
are applied for the smooth function Mhu: 



\MhU-IilMhu\\^rk,p(^Q^)<s^'ze{Qi) ^ ^ \MhU -IliMhu\1yt:, 

[cos Ok f J 



Lemma 13.61 is applied to combine terms in the summation and Lemma 13.31 is used to return the 
estimate to the correct norm. Letting 9i = maxTnQi^0^T, 

< size(Q.)'? ((^^) 

~ ((cos'^oO ■ 

Inequality ()30p follows from Proposition 13.51 which guarantees that h^^size{Qi) < 1. Finally sum- 
mation over all cubes and again utilizing Proposition 13.51 finishes the result: 

he 



I I IP 



E 

js.t.QjnQ,7^0 



(cos^O^^ 



The final inequality results from the fact that there are fewer than 4^ cubes intersecting Qi and for 
any such cube Qj, hq. < 2/iq-. □ 

Remark 3.G. Due to Remark 3.E, the k = case of Theorem 13.71 does not depend upon (I16p . 
This means that error estimates can be established for Lipschitz sizing functions (which, recalling 
Remark 3.C, result from adaptive refinement in many common settings). 

Remark 3.H. Similarly to Theorem 12. 4^ the extension domain property is also important to 
Theorem 13.71 The usage is somewhat implicit: the continuous extension property is used to define 
u on Qi <f_ Vt. 
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Figure 9: Requiring a globally defined sizing 
function can lead to restrictions associated with 
non-local portions of the boundary; e.g., a small 
mesh in region A will force a reasonably small 
mesh also in region B. 



4 Final Remarks 

To extend the maximum angle condition and its generalizations to less regular data Section [3] es- 
tablishes ply-independent error estimates for a linear average interpolant. Broadly, this serves to 
demonstrate that requirements on shape regularity and data regularity can be weakened simulta- 
neously without the inherent trade-off seen in previous analyses. Specifically, these results extend 
optimal error estimates under the maximum angle condition to several situations for which it was 
previously not shown: linear interpolation of functions in W^'^{Q) and (when 1 < d/p) W'^'P{Q). 
The latter result is most important in three dimensions where the maximum angle condition does 
not hold for the Lagrange interpolant when p <2. 

We have not addressed the question of boundary conditions. For error estimates in the L^- 
norm, the function can be explicitly adjusted to match a homogeneous boundary condition. In [11] 
averaging regions near the boundary are shifted outside the domain where an explicit extension 
by zero is used. When only considering the more straightforward scalar functions, an explicit 
cutoff can also accomplish the task \3^. However, these manipulations corrupt error estimates 
of the derivatives and fully generalizing (j28p to preserve boundary data requires an alternative 
approach. Ideally a mollification procedure analogous to the Scott and Zhang interpolant [38] can 
be constructed which averages over a lower dimensional manifold so that only boundary data is 
taken into account. Nai've approaches in this direction have proven unsuccessful. 

The extension domain property has been used somewhat implicitly throughout the arguments 
presented posing some limitations on the set of admissible meshes. Since the sizing function is 
defined globally, mesh size in non-local portions of the domain can impact the maximum mesh 
size; see Figure [H These limitations are often ignored in the analysis as asymptotically (as h 
becomes small) the impact disappears since the domain is fixed. Other mollification-based average 
interpolants face similar limitations. For example, the construction of Christiansen and Winther 
explicitly construct a small extension region outside the domain boundary and mesh elements are 
assumed to be smaller than the size of this extension region. Ideally the domain could be embedded 
in a manifold in which the non-local portions are far away (or more technically the extension 
operator has a much smaller constant). While ad hoc methods for doing this are apparent for 
simple examples, a complete theory for general domains is not. 

In some cases, a higher-order variant of Theorem 13.71 may be necessary. Recall that for the 
expected convergence rate in the VF^'^-norm, Jamet's theorem requires k + l — m > ^, where fc + 1 
is the order of interpolation used and d is the spatial dimension. Thus average interpolation is 
needed when d is large, p is near 1, or m is near k, i.e., estimates on higher derivatives are desired. 
Extending Theorem 13.71 to cover these cases requires a higher-order mollification procedure. While 
this mollifier can be defined, precise technical requirements under which the interpolation error 
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estimate holds are still unknown. 
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A Coplanarity vs. the Maximum Angle Condition 

Proof. (Proposition II. ip For this proof we omit the subscript K on 9k and without ambiguity. 
First, we will bound ip from above by a function of 9 to establish that if 9 is bounded away from 7r/2, 
then Ip is bounded away from tt. This will be shown by direct construction: for a tetrahedron with 
maximum angle ip, we will construct a particular vector ^ which ensures that 9 is sufficiently large. 
Construction is divided into two cases. First we suppose that ip is realized by a dihedral angle of the 
tetrahedron. Considering this largest dihedral angle, all vertices of the tetrahedron are contained 
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in two planes that meet at angles of "0 and vr — ■0. Selecting a vector that points "away" from these 
planes demonstrates 

> (31) 

Next we suppose that ^ip is realized by an angle in one of the faces of the tetrahedron. Let vi be 
the vertex at which the angle ijj occurs, let V2 and V3 denote the other two vertices of that triangle 
and let V4 denote the last vertex; see Fig. [TOj Moreover, without loss of generality we can assume 
that |v2 — V4I < |v3 — V4I, i.e., V4 is nearer to V2 than V3. Let ^ be the unit normal vector to the 
triangle with vertices vi, V2 and V4. To bound 9, first observe that $, is orthogonal to three of the 
six lines containing edges of the tetrahedron. Two other edges belong to the triangle with vertices 
vi, V2, V3, i.e., the triangle with a largest angle V'- This large angle ensures that for any unit vector 
Uj in the direction of either of these two lines, ^ • Uj > a — 7r/2. Lastly, we consider unit vectors in 
the direction of the line through V3 and V4. The triangle inequality implies 



^ • (V3 - V4) 


< 


i- (V3 - V2) 


+ 


1 • (V2 - V4) 




^ • (V3 - V2) 



where the second term in the sum is zero by construction of ^ orthogonal to the a plane containing 



V2 and V4. Recalling our earlier argument, 



IV3 - V2, 



< |v3 — V2I COS (V' — vr/2). Thus we can 



estimate the angle between ^ and the line containing V3 and V4. 



^ • (V3 - V4) 



< 



IV3 - V4I 



|v3 — V2I cos {iIj — vr/2) 
|V3 - V4I 



< 2 cos (V'-7r/2). 



(32) 



The final inequality above follows from our definition of V3: |v3 — V2I < |v3 — V4I + |v4 — V2I < 
2 |v3 — V4I. We thus conclude that 



> min^ • Uj > arccos (2 cos {ip — vr/2)) . 



(33) 



Considering the two expressions ()3ip and ()33p . we see that if 9 is bounded away from 7r/2, then ^ip 
is bounded away from vr. 

Next we establish a reverse inequality. The crux of the argument is contained in an inequality 
from [2j. Defining m := sin^^^^^, Acosta et al. assert that there are unit vectors Ui, U2, U3 in the 
directions of the edges of the tetrahedron such that 

|(ui XU2)-U3| >m3. [21 p. 147] (34) 

Letting ^ be an arbitrary unit vector, we consider the decomposition ^ = a(ui x U2) + 6ui2 where 
U12 is a unit vector in the plane spanned by ui and U2. Since these are orthgonal subspaces, 
a2 + &2 = 1 and thus a > I/V2 or 6 > 1/^/2. If a > I/V2, then 

n3 



1^ • U3I > a |(ui X U2) • U3I > 



7f 



(35) 



Otherwise b > l/\/2- Letting 7 denote the non-acute angle between the lines containing ui and U2, 
we observe that (|34p implies that sin (vr — 7) > or 7 < vr — arcsin(m^). Also the angle between 
U12 and the line in the direction of either ui or U2 is no larger than 7/2. Thus, 

,3^ 



max 1^ • Uj 

je{i,2} 



> 6 cos 



1 / vr 

> cos 
" V2 



arcsm m 



(36) 
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Combining ([35]) and ([36]) . 



. / f 1 /vr — arcsin(m^) 

& = max mm arccos £ • u,- < arccos max < —= , —= cos 

l€l=i * V l\/2'^/2 V 2 

Despite the complex form, this is the desired inequahty. Specifically, if ip is bounded away from vr, 
then m is bounded away from zero. Since all the functions we consider are continuous, it follows 
that the argument of arccos is also bounded away from zero. Thus 9 is bounded away from 7r/2. □ 




B Mollification 

Proofs of several propositions in Section 12.11 are provided here. 

Proof. (Proposition 12. 1[ ) Define by | + ^ = 1 and let w G Lp'(M°') be an arbitrary function. If 
u G C^{R'^) then 



/ \{u{x)-Mhu{x))w{x)\dyi 

(x) - / u(yL-y)ph{y)dy] w(x) 



dx 



/ (/ {u{x.) -u(yL-y))phiy)dy] wix) 



dx. 



The final equality results from the fundamental theorem of calculus. In the support of moUifier p^, 
\y\ < h so \yph{y)\ < hph{y). Then Holder's inequality can be applied: 



\{u{-x.)-Mhu{x.))w{-K)\ dx 

/ / / |-y • Vn(x - sy)p/i(y)tt;(x)| dxdy ds 

Jo JR'i JR'i 

|tt;(x)|^ dx ) dy ds 



dx 



< 



< / / hph{y){ I |V^z(x-sy)|Pdx 
Jo Jr"^ 

h ||Vti||^p 



By selecting tt; = {u — M-^uY ^, the left-hand side of the above estimate becomes \ \u — A^/iu||^p 

Since p'{p- 1) = p, \\w\\Lp'{R<i) = 11^ ~ -^/i'"llLP(]Rd)- Dividing both sides by | lip'(Kd) gives the 
desired inequality. Density of smooth functions in is used to extend this result to the full 

space. □ 
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Proof. (Proposition I2.2p Let u{x) := u{x)xq,S^) where XQh the characteristic function of Qh- 
Applying Young's inequaUty to M.hU gives 

Since MhU = MhU on Q, \\Mhu\\i^p^Q) = \\Mhu\\LP(Q) < \\-Mhu\\LP{Rdy Then recognizing that 
II^IIl9(ir'*) ~ computing the norm of ph yields the desired estimate. □ 
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Figure 10: In the proof of Proposition II. H when 
bounding ■0 in terms of 9, a specific $, can be con- 
structed explicitly when -0 is realized by a large di- 
hedral angle (left) or a large angle in a face (right). 
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